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We analyze the mean-field eigenstates of the atom-molecule Bose-Einstein condensate (AMBEC) 
under the assumption that the background (elastic) scattering length of the atoms can be ignored. 
It is shown that the relevant eigenstates are localized in the space of the condensate parameters: The 
eigenstate has a different character in different regions of the parameter space, and at the interface of 
two local eigenstates the properties of the system may change nonanalytically. Using both analytical 
and numerical techniques, we find the approximate boundaries of the local eigenstates and identify 
the types of parametric excitations that occur when an eigenstate is forced outside of its region of 
validity by a parameter sweep. We contrast the properties of the mean-field parametric excitations 
found in AMBEC with the experimentally observed excitations of the BEC. 
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I. INTRODUCTION 

The dynamics of the dilute, trapped Bose-Einstein 
condensates (BEC) of a single atomic species has been 
successfully modeled using the Gross-Pitaevskii equation 
(GPE) [1, 2] [toss these reference], 

= (-^V 2 + V (r,t))$M) 
+ ^^|$(r,i)| 2 0>(r,t). 

Here $ is the mean-field amplitude of the atomic con- 
densate, normalized to unity, N is the number of atoms 
in the condensate, m is the mass of an atom, Vq is the 
trapping potential, and a is the atom-atom s-wave scat- 
tering length. The properties of the mean-field ground 
state <E> of Eq. (1) depend on the sign of a. If a > the 
ground state $ is a single, well behaved function. How- 
ever, when a < 0, there may be two solutions: the ground 
state, which is the collapsed state of the condensate rep- 
resented by a (5-function, and a metastable state, which 
acts as a continuation of the positive-a ground state to 
the negative values of a. The metastable eigenstate exists 
only for limited values of the parameters N a/i ~ —0.67, 
where £ is the characteristic length of the harmonic trap- 
ping potential Vq. If N a/l is driven beyond this limit by 
increasing, say, the atom number in the condensate, the 
amplitude <E> collapses [3-7] . It has been argued that the 
collapse in the BEC takes a finite time [8] and leads to a 
gas- liquid phase transition [9-11]. 

Controlling the collapse is but one motivation for the 
current flurry of effort toward controlling the scatter- 
ing length [6, 12-14]. Generally, atom-atom interactions 
can be manipulated experimentally by adjusting the cou- 
pling of a pair of atoms to the corresponding molecular 
dimer with a static magnetic field (Feshbach resonance) 
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or with electromagnetic radiation (photoassociation) [15- 
19]. [add Stwalley's suggestion]. These schemes may be 
modeled with a quantum field theory of coupled atomic 
and molecular fields. A formally identical field theory can 
be devised for elastic resonant scattering [20], with the 
difference that the scattering resonances takes the role of 
the molecular dimer. All of these approaches lead to the 
same mean-field theory of the atom-molecule BEC (AM- 
BEC) in which there are three parameters: K, coupling 
strength between pairs of atoms and molecules; S, detun- 
ing of the molecular state from the edge of the two-atom 
continuum; and a^g, background atom- atom scattering 
length. 

In this paper we examine the eigenstates of the mean- 
field AMBEC under the simplifying assumption that 
o-bg = 0. The focus is on two eigenstates: the thermody- 
namically stable ground state, and a "twin" state that, 
we believe, is often more relevant than the ground state 
in experiments starting with an atomic condensate. Us- 
ing numerical and analytical techniques we show that the 
{5, K} parameter space is divided into regions, between 
which the nature of the eigenstate abruptly changes. We 
associate the existence of such "local" eigenstates within 
a "global" eigenstate with the non-linearity of the mean- 
field theory of the coupling between the atoms and the 
molecules. 

When the system is created in a local eigenstate and 
is subsequently pushed by a time-dependent variation of 
the parameters into a region in which the initial local 
eigenstate no longer exists, parametric excitations ensue. 
We will characterize and classify the excitation analyti- 
cally and numerically. 

Finally, we will point out intriguing similarities with 
experiments. Caution is due in such comparisons, as the 
motion of the atoms and the molecules in the trap is 
significant in our calculations and not necessarily so in 
all experiments. Conversely, though, in genuine trap ex- 
periments heretofore unrecognized possibilities open up. 
In particular, with a proper variation of the parameters 
the collapse of a condensate with a negative (effective) 
scattering length may be achieved continuously and in a 
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controllable fashion. 



III. EIGENSTATES OF ZERO-DIMENSIONAL 
AMBEC: GROUND STATE AND TWIN STATE 



II. MEAN FIELD MODEL OF SIMPLE AMBEC 

We use the mean fields ip and ip to describe atomic and 
molecular condensates, respectively. The interaction en- 
ergy produced by the coupling between the condensates 
reads 



(2) 



In the theory of Feshbach resonances, K 2 is proportional 
to the width of the resonance [21-24]. Conversely, in 
photoassociation, K 2 is a controllable parameter propor- 
tional to laser intensity [18]. 

Our assumptions of the simple AMBEC lead to the fol- 
lowing set of scaled coupled Gross-Pitaevskii equations, 



i^- = H a ip - Kip*ip, 
&- = (H m -5)iP-K^ 2 . 

OT 



(3a) 
(3b) 



The detuning 8 measures the difference in energy of an 
atomic pair and the discrete molecular state. Here we 
choose the sign of <5 so that a positive (negative) detun- 
ings produces a positive (negative) effective scattering 
length. The Hamiltonians H a and H m correspond to the 
center-of-mass energies for atoms and molecules, respec- 
tively, 

H a = -iv 2 + i]T<^, (4a) 

Hrn = -V+^<^ 2 . (4b) 



For simplicity, we assume that the trap is isotropic with 
cj = 1. Lastly, we scale the fields so that the normaliza- 
tion reads 



(<p\<p) + (V#> = 1. 



(5) 



The expression for the conserved total energy of the sys- 
tem is then 

E=(<p\H a \<p) + ±{il>\H m -6\il>)+Ei. (6) 

There are numerous practical issues with the model 
embodied in Eqs. (3), many of which have been discussed 
elsewhere, e.g., in Ref. [18]. Here we focus on the AM- 
BEC eigenstates, either global (stretching over the whole 
parameter space) or local (limited to certain regions of 
the parameter space), an important issue of principle as- 
sociated with variations of the AMBEC parameters 8 and 
K. For the most part, we ignore the complications of real 
experiments. 



The system of CPE (3) of a simple AMBEC is non- 
linear in the amplitudes ip and ip. For that reason we start 
our discussion of the AMBEC spectrum with an analysis 
of its simpler, zero-dimensional, version, and recite some 
elementary results [15, 18, 25]. 

A zero-dimensional AMBEC model is obtained by re- 
placing the Hamiltonians H a and H m in (3) with con- 
stants, which are subsequently absorbed in 8. The time 
evolution of the atomic (a) and molecular (/3) amplitudes 
is then governed by a system of ordinary differential equa- 
tions (ODE), 



ia — —Ka*j3, 
i/3 = -8/3 -Ka 2 . 



(7a) 
(7b) 



The stationary solutions arc sought in the form a = 
x ■ exp(— i [i t) and fi — y ■ exp(— i 2 fit). Additionally, the 
normalization (5) transforms into x 2 + y 2 = 1 . 

We find three possible solutions: the ground state, the 
twin state and the all-molecule state. The all-molecule 
state is a trivial solution with the amplitude y = 1 and 
the eigenenergy \i = —8/2, where 8 — 8/K. The ground 
state has the molecular amplitude 



yas(S) = 



and the frequency is 




for 6_ < 2, 
for 5 > 2, 



(8) 
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for 6 < 2, 
for 6 > 2. 



(9) 



The twin state is specified by 



Vtw(S) 



1, 



6 ■ 



for 8 < -2, 
for 8 > 2, 



with 



I^Tw(d) 



\A 2 + 12 

6 - 



for 8 < -2, 
for 8 > 2. 



(10) 



(11) 



In terms of the fraction of atoms, the ground state (8) 
behaves as follows. For large negative detunings, 8 <C — 1, 
the ground state consists almost entirely of atoms. The 
fraction of molecules increases as we approach 8 C = 2 
from below. At the non-analytic point 8 = 8 C the ground 
state merges with the all-molecule state (y = 1), but that 
merge is not smooth. The twin state (10) behaves sim- 
ilarly, as the ground and twin states transform between 
each other with the operations 8 — > — 8 and y — > — y. For 
large positive detunings, 8 S> K, the fraction of atoms 
is almost unity and it decreases as we approach 8 = — 2 
from above. At 8 — —S c the twin state merges with the 
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all-molecule state (y = — 1). However, as was the case 
with the ground state, the transition to an all-molecule 
state is continuous but not smooth. 

We next examine the stability of the stationary solu- 
tions of Eq. (7) by finding their frequencies of small oscil- 
lations. The most convenient variables are the molecular 
amplitude y and the phase difference between the square 
of the atomic amplitude and the molecular amplitude, 
9 = 2arg(<p) — arg(?/>). In these new variables the origi- 
nal system of ODE's (7) reads 



y 



(l-y 2 )sin0, 



= S + ( 3y ) cos ( 



(12a) 
(12b) 



Its nontrivial (x ^ 0) fixed points are constructed from 
the solutions (8) and (10) so that the amplitude y is 
positive. The ground state is then y = yosi (8), with 
phase 9 = 0. The twin state has the phase 9 — it and 
consequently y = —yr, (10). Both stationary solutions 
have the same frequency of small oscillations derived from 



n 2 = cos 2 (1 



•3), 



(13) 



which always gives a real number. 

For the all-molecule solution the analysis is somewhat 
more complicated. For \S\ < 2 the all-molecule state is 
unstable [25], while for \S\ > 2 it is stable. To get a feel 
for how the stably the all- molecule state behaves, let us 
consider the large 6 — S/K limit and write the molecular 
amplitude as y = 1 — Sy. In that limit the approximate 
solution for phase 9 is 9 s=a St. Solving (12a) for 5y(t) 
yields 



Sy(t) w Sy(0) ■ exp 



V2 



(1 - cos(St)) 



(14) 



That is, for \S\ > 2 the all-molecule solution is stable, 
but not an attractor: the oscillations of amplitudes in 
the vicinity of the all-molecule state are persistent and 
do not decay with time. 

To conclude, the zero-dimensional model (7) has three 
stationary states: the ground state, the twin state and 
the all-molecule state. The ground state (twin state) has 
a non-analytic point at 5 — +2 (—2) where it merges with 
the all-molecule state. A detuning sweep across either 
of the non-analytic points creates parametric excitations 
of the amplitudes, the details of which depend on the 
numerical details (sweep rate, initial conditions). Com- 
monly, these excitations manifest themselves in small, 
persistent, non-linear oscillations near the all-molecule 
state. 



IV. GROUND STATE AND TWIN STATE IN 
SIMPLE AMBEC - NUMERICAL STUDY 

The simple AMBEC (3) inherits much of the behav- 
ior of the zero-dimensional version, but also adds some 



unexpected twists. Before going into the details, we ex- 
pand on some terminology we have already brought up 
in passing. 

The zero-dimensional model has stationary states, or 
eigenstates. A notable one is the ground state, the eigen- 
state with the lowest chemical potential. There are no 
general criteria for the existence of eigenstates for non- 
linear differential operators, but it appears from our nu- 
merical calculations that the simple AMBEC, too, has 
a unique ground state for all parameters S and K. As 
the ground state can be identified, seemingly unambigu- 
ously, in the entire parameter space, we call it a global 
eigenstate. 

In the zero-dimensional model the nature of the ground 
state is different depending on the parameters. For in- 
stance, if S > 2, the ground state is all molecules. A 
similar behavior is found in the ground state of the AM- 
BEC. Locally, in different regions of the parameter space, 
the ground state may be (nearly) all molecules, (nearly) 
all atoms, or a mixture of atoms and molecules. As in the 
zero-dimensional case, the nature of the ground state may 
change nonanalytically when the AMBEC parameters are 
varied. We say that at such a point of nonanalyticity the 
global ground state of the AMBEC switches from one 
local (in the {S, K} space) eigenstate to another. We 
characterize local eigenstates as all-atom, all-molecule, 
or mixed atom-molecule. On occasion, this distinction 
is qualitative only. For instance, the all-atom state may 
not be all atoms, but the fraction of atoms is much larger 
than in the adjacent mixed atom- molecule state. 

However, the experiments do not always deal with the 
ground state. In a Feshbach resonant system, what seem- 
ingly is an atomic BEC may be prepared with parameters 
such that the ground state of the AMBEC would be all 
molecules. We introduce the twin state analogous of the 
twin state in the zero-dimensional case to model this type 
of a situation. We construct the twin state numerically 
by starting with an all-atom ground state of the nonin- 
teracting system (K = 0), and then vary K adiabatically. 
This procedure appears to produce an eigenstate of the 
AMBEC system for all 5 > and K . It appears that for 
certain K > and 5 < the twin state is unstable, and 
in practice it is not possible to construct it numerically 
by integrating the AMBEC system in time. Nevertheless, 
we surmise that we have here another global eigenstate 
of the AMBEC. Just like the ground state, we expect the 
twin state to be split into three local eigenstates,: all- 
atom, all-molecule, and mixed atom-molecule eigenstate. 

Except for stationary states, we will also investigate 
what happens when the parameters 8 and K are varied in 
time. The general observation is that, where the system 
switches from one local eigenstate to another, various 
type parametric excitations set in. 

The present Section IV reports on numerical studies 
about the breakdown of the ground state and the twin 
state into local eigenstates, and describes the parametric 
excitations that ensue when a parameter is swept across 
a border between local eigenstates. 
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A. Ground State 



Here we calculate the ground state numerically. While 
the work in three dimensions is feasible in principle, the 
computations are restricted to spherically symmetric sit- 
uations for simplicity. We employ the DS- method [26]. 
The iteration that is repeated until convergence consists 
of integration in complex time followed by the renormal- 
ization, see Appendix A. For a point (6, K) in parameter 
space we find the solutions ip and ip for the amplitudes. 
Next we determine the atomic fraction, which is the same 
as the norm of atomic amplitude, N a — (<p\ip), and the 
half-size J?i/2, Eq. (A4), of the atomic distribution. The 
calculation of the ground state amplitudes is done over 
an integer mesh S = —30... 30, and K = 0.1,1... 30; 
where we use K = 0.1 as an approximation of the limit 
K -» 0. 

The fraction of the atoms and the half-size in the 
ground state are shown in Figs. 1 and 2, respectively, 
for different points in the {8, K} parameter space. We 
observe that the ground state consists of three local eigen- 
states: all-atom state, mixed or atom-molecule state, and 
all-molecule state. We refer to the boundaries of partic- 
ular eigenstate as a fractures, and label them as follows. 
The all-atom state is bounded by the fracture Jq S from 
the mixed or atom-molecule state and similarly, the all- 
molecule eigenstate is bounded by the fracture Jq S from 
the mixed or atom-molecule state. 

We next demonstrate how a parameter sweep of the 
eigenstate across the fracture creates a parametric ex- 
citation. In the first example we perform an intensity 
sweep, with K = . . . 18 and 5 = —30. The behavior of 
the atomic fraction and the half-size is shown in Fig. 3. 
As surmised, at the fracture Jq S the local all-atom eigen- 
state disappears causing a parametric excitation of both 
amplitudes. The extent of the excitation is determined 
by the relative position of the mixed or atom-molecule 
eigenstate (which becomes a new ground state) and the 
numerical details of the sweep (rate, initial conditions) . 

In the second example we perform a positive detuning 
sweep, with 5 — —30 ... 30 and K = 10. Fig. 4 shows 
the behavior of the norm and the half-size of the atomic 
distribution. For 5 < we find no evidence for para- 
metric excitations in the amplitudes, meaning that the 
all-atom state may continuously evolve into the mixed or 
atom-molecule state in some regions of parameter space. 
For S > 0, close to the fracture Jq S , we observe another 
parametric excitation. 

The parametric excitations at Jq S and Jq S differ. At 
Jq S the oscillations of the atomic quantities are erratic, 
whereas at /q S the oscillations are periodic. We return to 
these differences in Sec. V, when we analyze the solutions 
for the amplitudes analytically. 



B. Twin State 

In the zero-dimensional model, the twin state and the 
ground state are related via the transformation 5 — > 
—6,K — ► —K. Applying this transformation to the 
ground state, we surmise that the twin state contains 
the all-atom state for S > 0, K — > 0, and may include the 
all-molecule state. 

We test our assumptions using numerical methods. For 
parameters K > and S > 0, the twin state is prepared 
from the all-atom state, ground state of the atomic sys- 
tem in the given trap, by slowly varying the matrix ele- 
ment K from zero to the desired value at the fixed de- 
tuning 5. Following such an adiabatic preparation, the 
twin state is monitored in the standard fashion while one 
of the parameters, either the detuning S or the matrix 
element K, is varied. 

We limit the exposition of our results to two charac- 
teristic examples. In the first example, shown in Fig. 5, 
we fix the detuning at 5 = 100 and vary K in the range 
... 36. As can be seen, the fraction of the atoms stays 
close to 1 at all times, while the half-size of the atomic 
condensate increases slightly. We observe that there are 
no parametric excitations, and conclude that the all-atom 
state exists for 8 > and K > 0. This reaffirms adiabatic 
intensity sweep as a method for preparing the twin state 
for positive detunings, and demonstrates the absence of 
the analogue of the ground state fracture Jq S in the twin 
state. 

In the second example we take the twin state prepared 
at 5 = 50 and K = 10, and perform a negative detun- 
ing sweep 5 — 50 ... — 50. The behavior of the norm 
and the half-size of the atomic distribution are shown in 
Fig. 6. For 6 > the all- atom state has no parametric 
excitations, in accord with the first twin-state example. 
Around S w the fraction of the molecules starts to in- 
crease and the all-atom eigenstate becomes a mixed or 
atom-molecule eigenstate with an increasing fraction of 
molecules. At some negative S we find familiar signatures 
of a strong parametric excitation in the atomic quanti- 
ties, which suddenly start to oscillate chaotically. We 
refer to the set of points in the {5, K} parameter space 
where the parametric excitations set in as the fracture 
frw- 



V. ANALYTIC BOUNDARIES OF LOCAL 
EIGENSTATES 

Both the twin state and the ground state have as a lo- 
cal eigenstate either an all-molecule state (zero-fraction 
of atoms), or an all-atom state (negligible fraction of 
molecules). The boundaries of those two special eigen- 
states in the parameter space may be found using an- 
alytic methods. Our next task thus is to find analytic 
expressions for the boundaries of the local eigenstates. 
We also compare the results with direct numerical com- 
putations. 
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FIG. 2: Half-size R1/2, Eq. (A4), of the atomic distribution in the ground state as a function of the detuning 8 and the 
atom-molecule coupling K. 



A. All-atom state; fracture /q S 

In the limit K — > and » Jf the fraction of 
molecules becomes negligible, so it is possible to elimi- 
nate the molecules from the theory. Performed in the 
standard way, this procedure leads to a single-condensate 
GPE with a negative scattering length. Its solution is ci- 
ther a single collapsed state or a pair consisting of the col- 
lapsed state and the "mctastablc" (non-collapsed) state. 
The collapse boundary in this approach is the boundary 
of the metastable state in the parameter space. In the 
variational approach [3] the collapse boundary is given by 



a c w —0.67, which in terms of the AMBEC parameters 
reads w -0.67. 



We observe that there is a contradiction between 
the "no-collapse" theorem [27], according to which the 
AMBEC amplitudes may not collapse, and the "stan- 
dard" single-condensate approximation of the AMBEC, 
in which the collapse of the atomic amplitude is allowed. 
We resolve the contradiction by reexamining the effective 
scattering length expansion, and demonstrate that the 
usual expansion is incomplete. We recall that Eq. (3b) 
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FIG. 3: (Ground State) Norm N a and the half-size R1/2 of 
the atomic distribution for the intensity sweep K — ... 18 
with fixed 8 — —30. Integration is performed in 50,000 steps 
of size (A2). At f GS the atomic variables begin to oscillate. 
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FIG. 4: (Ground State) Norm N a and the half-size R1/2 of 
the atomic distribution for the detuning sweep S — —30 ... 30 
with fixed K = 10, done in 50,000 steps of size (A2). Fol- 
lowing Jq S the half-size of the atomic condensate begins to 
oscillate erratically. 
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FIG. 5: (Twin State) Norm N a and the half-size R 1/2 of the 
atomic distribution for the intensity sweep K = . . . 36 with 
fixed <5 = 100, done in 150,000 steps of size (A2). We see no 
evidence of parametric excitations. 



FIG. 6: (Twin State) Norm N a and the half-size R 1/2 of 
the atomic condensate for the negative detuning sweep 8 = 
50 ... - 20 with fixed K = 10, done in 100,000 steps of size 
(A2). The time is reversed (goes from right to left) so that the 
detuning axis has the same orientation as in Fig. 4. Following 
frw the molecules decay rapidly into atoms and both atomic 
and molecular distributions begin to oscillate erratically in 
half-size and fraction. 



can be formally solved for the field tp [27] , yielding 



!)l 



6-Hr, 



_K v^oo n / . d TT 

<5 2^j=0 IS Vdt n " 



IV 



(15) 



In the limit N m = (ip\tp) <C 1 and large \S\, the term 
iip can be neglected in (15) when com- 



at 



pared to H m ip ~ ^ip. 

After removing the partial time derivative in (15), we 
group the terms with j = and j = 1, then j = 2 and 
j = 3, and so forth. The new expansion is well behaved; 
if the amplitudes collapse for whatever reason, <p(r) — > 
<5(r), each term in the rearranged expansion grows to 
+00. The 6- function is thus a solution with an infinite 
positive energy and the collapse of the condensate is not 
energetically favorable, in accord with the "no-collapse" 
theorem. The correct leading-order approximation for 
H' a , which combines the terms j = (effective scattering) 
and j ' = 1 (formation of a coherent atom pair) , reads 

ff> = ^ + ^MV + §jV 2 iW, (16) 

with H a and H m from Eqs. (4). We observe that in this 
approximation a molecule (a new entity made of two 
atoms) is replaced by a quasi-particle (a coherent excited 
state of two atoms) . 

We examine the eigenstates of the Hamiltonian H' a us- 
ing the variational approximation [15, 28, 29]. The total 
energy is given by 

E a {fp) = (ip\H a \<p) + ^ (^V) + (<P 2 \H m \<p 2 ) • 



(17) 



Here the atomic eigenstate ip is approximated by a Gaus- 
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TABLE I: Analytical (K t h) and numerical (K num ) positions 
of f GS in the parameter space. The atomic fraction (N aa ) 
and half-size of the atomic distribution (-R1/2) are evaluated 
numerically on the opposite sides of fas- f° r the all-atom 
eigenstate at K — K num — 0.5 (aa), and for the mixed or 
atom molecule state at K = K num + 0.5 (mix). 



FIG. 7: Variational sizes of atomic distribution for the detun- 
ings 8 = -10 and 8 = -50, with K = ... 30. For 8 = -10, 
there is a single non-physical solution. For 8 = —50, there are 
two solutions, one of which is the ground state. 
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FIG. 8: Adiabatic evolution of the initial all-atom eigenstate 
at the fixed detuning 8 — —5 for K = . . . 15, in terms of the 
fraction of atoms and the half-size of the atomic distribution. 
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FIG. 9: The position of f GS in the parameter space, found 
using the variational approximation. 
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(18) 



The total energy becomes a function of the variational 
condensate size s. The variational eigenstates are the 
states of the form (18) for which the size s is a local 
minimum of the energy functional E a (s). That is, the 
variational sizes {si}i = i... are the real positive roots of 
dE/ds(si) = such that d 2 E/ds 2 (s l ) > 0. We use 
Mathematica© to find the solutions for the variational 
size of the atomic condensate, and keep only those solu- 
tions which do not exist in the limit K — > oo for a given 
5 < 0. Our choice of a physical all-atom eigenstate is jus- 
tified considering the behavior of the fraction of atoms in 
the ground state shown in Fig. 1. In the limit K — > oo the 
ground state is not an all-atom state. In this approxima- 
tion, for given 8 < 0, the fracture /q S occurs at the value 
of K for which the all-atom ground state disappears. 

To illustrate this point, let us consider the variational 
size for two fixed detunings, 6 = — 10 and 5 = —50, and 
K = 0...30, shown in Fig. 7. For <5 = —10, there is 
a single solution which exists for all K, thus it is not 
acceptable as a solution. For 5 = —50, there are two 
solutions. The lower one exists in the limit K — > oo, so 
it is not acceptable. For the true ground state we then 
take the other solution that exists in the limit K — > 0. At 
the critical detuning 8 C ~ —17.8, the number of solutions 
switches from one (one physical and one not) to zero (one 
non-physical) . 

The fact that /q S only partially separates the all-atom 
state from a mixed or atom-molecule state is illustrated 
in Fig. 8. We observe a smooth decrease of the fraction 
and the half-size of the atomic distribution by increasing 
K from to 15 at a fixed S — —5. 

We plot the analytical values for Jq S based on our vari- 
ational approximation in Fig. 9. The tip of the fracture 
occurs for 5 w —17.8 and K w 13.5. 

To determine the position of the fracture numerically 
for a fixed S, we have covered the parameter K in in- 



8 




FIG. 10: Outline ol the variational solutions ol A = and 
X TW = 0. In regions A > 0, there is no solution for the vari- 
ational size s. In regions where A < there are two solutions 
for the variational size s. Fractures frw and f GS appear as 
the boundaries between the regions. 
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til/2 


-1.5 













2.6 


2.5 ± 0.5 


0.13 0.76 


1 


4.1 


4.5 ± 0.5 


0.18 0.76 


5 


9.1 


9.5 ± 0.5 


0.25 0.52 


10 


14.5 


14.5 ± 0.5 


0.42 0.28 


20 


24.2 


24.5 ± 0.5 


0.55 0.09 


30 


33.4 


33.5 ± 0.5 


0.66 0.03 



TABLE II: Analytically (Kth) and numerically {K nu m) cal- 
culated fracture fas m t ne K} parameter space. The frac- 
tion of atoms and the size of the atomic distribution are given 
for the mixed or atom-molecule state (superscript mix) at 
K = Knum + 0.5. 



crements of 1, and assigned each K to one of the local 
eigenstates. In this manner we find the position of the 
fracture line to within ±0.5. In Table I we compare the 
position of the fracture found directly numerically, and 
analytically using the variational method. The agree- 
ment is very good. 

In Table I we also report the fraction of atoms and 
the half-size of the atomic distribution for the two adja- 
cent values of K that we have assigned to different lo- 
cal eigenstates. For instance, the half-size of the atomic 
distribution on the mixed or atom-molecule side of the 
fracture decreases when moving away from the origin. 
This change in size goes from i?™^ = 0.33 at S = —20 
to = 0.05 at 5 = -100, and more. While in AM- 

BEC there is no formal collapse of the amplitudes, for 
large negative values of 5 the substantial decrease in the 
half-size across Jq S can convincingly mimic a collapse. 



FIG. 11: Regions of A 2 , Eq. (28), in the parameter space. 
A 2 is positive only in Region III, bounded by fasi for all 
variational atomic sizes s. In all other regions there exist 
unstable variational solutions for s, for which A 2 < 0. 

B. All-molecule state; fractures f GS and frw 

The all-molecule state is the trivial global eigenstate of 
the AMBEC (3). Numerics shows that the all- molecule 
state is part of the ground state, and strongly influences 
the twin state. We analyze the all-molecule state from 
the perspective of the existence and stability of infinites- 
imal atomic configurations. 

The amplitude of the all-molecule state is simply the 
ground state wave function of the molecular harmonic 
potential, 

^)=(^) 3/ V— 2 , (19) 

where m is the mass of the atom and uj is the trapping fre- 
quency. The norm of the molecular amplitude is unity, in 
accord with the all-molecule hypothesis, and the eigenen- 
ergy is 

* = 5 (§-'). (20) 

We now analyze the presence of atoms in a perturbative 
fashion. 



1. Existence of variational solutions for the atoms 

In the all-molecule approximation an infinitesimal frac- 
tion of atoms is roughly described by Eq. (3a), where 

W = H a <p-Kip*i>. (21) 

Here, for the atomic amplitude ip we use the variational 
expression from Eq. (18), with the size s as the varia- 
tional parameter. In (21), the eigenenergy (i is given by 
Eq. (20). Let us introduce a new quantity \ GS (s), 

X GS (s) = (ip\H a \ip)- f x-K(i P 2 \^), (22) 
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$th K t h 




-1.5 
-3.54 5 
-4.98 10 
-7.07 20 
-8.66 30 
-11.11 50 
-15.42 100 


-3.0 0.21 
-3.9 0.32 



TABLE III: Position of frw in parameter space in the varia- 
tional approximation. The position of the instability observed 
in the numerical simulations for a given K is shown in col- 
umn <5 ln . In column N™ is the calculated fraction of the 
atoms prior to parametric excitations. 



in terms of which the Eq. (21) becomes 
X GS (s) = 0. 



(23) 



In the expressions above the actual fraction of atoms can- 
cels out, so the norm of ip can be assumed to be 1. Anal- 
ysis of the real positive roots of Eq. (23) shows that there 
are two regions in parameter space, a mixed or atom- 
molecule region (A GS < 0), in which there are two vari- 
ational solutions for s, and the all-molecule region, in 
which there are no such solutions (X GS > 0). The frac- 
ture Jq S is the boundary between the two. This is shown 
in Fig. 10. 

This approach is easily modified to account for the all- 
molecule part of the twin state. From the zero-mode 
case, Eq. (10), it follows that in the twin state the phase 
shift between the molecular amplitude tjj and the matrix 
element K is —1. Allowing K to be negative (initially 
we said that K is non-negative) in Eq. (22) is a way to 
analyze the all-molecule state of the twin state. Again, 
we introduce a new quantity, X TW , 



A 



TW 



{^\H a \ v )-^ + K{^). 



(24) 



The variational size of the atoms s is the real positive 
solution of the equation 



X TW (s) = 0. 



(25) 



We perform the analysis of the real positive roots of 
Eq. (25), and find a behavior similar to the ground state. 
In the parameter space of the twin state there are two 
regions, the all-molecule region, in which there are two 

variational solutions for s (X TW < 0), and the mixed or 
atom-molecule region, in which there are no such solu- 
tions (X TW > 0). This classification is motivated by the 
numerical results, and is the opposite from the classifica- 
tion for the ground state given above. Nonetheless, the 
fracture f TW is the boundary between the two regions, 
as was the case in the ground state. Both fractures and 
the outline of the properties of the variational solutions 
are shown in Fig. 10. 

In Table II we compare the positions of the fracture 
Jqs found numerically (K num ) and analytically (Kth), 



for various detunings <5. Additionally, the table con- 
tains the fraction (N™ x ) and the half-size {Rffi) of 
the atomic distribution on the mixed or atom-molecule 
side of the fracture. Table II indicates that at /q S the 
two different local eigenstates, the all-molecule state and 
the mixed or atom-molecule state have the same (lowest) 
eigenenergy, but do not merge in terms of their properties 
(size, fraction, etc.). In Table III we compare the posi- 
tion of Jtw m the parameter space calculated using the 
properties of the variational solutions, Eq. (25), and the 
results from numerical simulations regarding the onset of 
parametric excitations. Numerical results are limited to 
two detunings, where for K = 10, 20 we find the detuning 
(S m ) and the fraction of atoms (A^ ra ) at the onset of the 
instability. 

We note that the analytical results for the all-molecule 
state differ from the (exact) numerical results in an im- 
portant aspect. Figs. 4 and 6 show the parametric ex- 
citations of the atomic amplitudes in the ground state 
and the twin state, respectively. In the ground state, 
a small atomic fraction survives the parametric excita- 
tion following Jq S , even though there are no available 
variational atomic configurations (X GS > 0). Conversely, 
the atomic fraction grows exponentially in what is sup- 
posed to be the all-molecule region of the twin state, 
although there are available variational atomic configu- 
rations (X TW < 0). We resolve this peculiar behavior by 
doing the stability analysis of the all-molecule state in 
the parameter space 



2. Stability of the all-molecule state 

The starting point for the stability analysis of the all- 
molecule state is the equation governing the behavior of 
the atoms, Eq. (3a). Let us write its solution in the form 



(/?(r) = (a(r) +ib(r)) ip exp(— i/ir), 



(26) 



where ip is the variational atomic amplitude, Eq. (18), 
and /i is the eigenenergy of the all-molecule state, 
Eq. (20). Here we assume that the time dependent real 
functions a(r) and b(r) are infinitesimals. We use this 
Ansatz in Eq. (21), multiply both sides with <p* and inte- 
grate in space. In this approximation Eq. (3a) is reduced 
to a system of ODE for a pair of real functions a and 6, 



b = -(^\H a \p)-^-K(^\4,)) a. 
We find the characteristic frequency of (27) to be 



(27) 



A 2 (s) = (MH a \v) - - K%^)\ 2 (28) 

As long as A 2 (s) > 0,Vs > 0, a small perturbation of the 
all-molecule state oscillates harmonically; and conversely, 
if 3s' > : A 2 (s') < 0, a deviation from the variational 
solution begins to grow exponentially. 
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The characteristic frequency is a product of two terms, 
A GS and X TW . The outline of their variational solu- 
tions in the parameter space has been given above and 
in Fig. 10. However, to have both fractures in the same 
half of the parameter space (K > 0), the fracture Jtw 
needs to undergo a reflection with respect to the K = 
axis (K -> -K). This told, in regard to Eq. (28) the 
parameter space is divided into three regions (I, II, III) 
separated by the fractures Jq S and frwi as shown in 
Fig. 11. In Region I, A 2 (s) = has 4 real positive solu- 
tions, so there are unstable atomic configurations. This 
region is bounded by Jtw- I n Region II, bounded by 
frw and fcs> ^ 2 ( s ) = has two real positive solutions, 
so there are unstable atomic configurations as well. It is 
only in Region III that A 2 (s) > 0, Vs > 0. In that part 
of the parameter space the all-molecule state is stable. 
Region I is what we call the all-molecule region of the 
twin state, while Region III is the all-molecule region of 
the ground state. 

We are now able to qualitatively describe the effect 
of an adiabatic sweep of the detuning on the ground 
and twin states. In a positive detuning sweep of the 
ground state, Fig. 4, the fraction of molecules increases 
while approaching the fracture Jq S . Somewhere in the 
vicinity of Jq S the mixed or atom-molecule state disap- 
pears and the only available close eigenstate is the all- 
molecule state. Disappearance of one eigenstate sets on 
the parametric excitations. The excitations manifests in 
non-periodic oscillations which are stable, because, as the 
stability analysis shows, the all- molecule state is stable 
in this region. However, because the atoms do not have 
an available variational size s, the shape of the atomic 
amplitude varies erratically in time. These variations are 
evident in oscillations of the half-size of the atomic dis- 
tribution, Fig. 4. 

The dynamics of the twin state under the negative 
detuning sweep follows a similar path prior to frwi as 
shown in Fig. 6. The fraction of molecules starts increas- 
ing toward frw, and somewhere close to the fracture the 
mixed or atom-molecule state disappears. At that point 
the only nearby eigenstate is the all-molecule eigenstate. 
However, because the all-molecule state is not stable in 
this part of the parameter space, the molecules start to 
decay into atoms. In the numerical simulations we did 
not find evidence for any other stable eigenstate for the 
system could settle in. The final result is chaotic os- 
cillations of the AMBEC amplitudes in both shape and 
fraction. 



VI. ABOUT THE EIGENSTATES 

In the previous section we calculated the boundaries 
of the local eigenstates by making use of their properties 
in the parameter space. We now reorganize these results 
into a description of the "global" eigenstates, the twin 
state and the ground state, to which the local eigenstates 
belong. 



A. Ground state 

The results of our discussion of the fractures in the 
ground state are compiled in Fig. 12, which shows the 
outline of the ground state in the parameter space. The 
ground state consists of two local eigenstates, the mixed 
or atom-molecule state and the all-molecule state, sepa- 
rated completely by fg S - Further, the mixed or atom- 
molecule eigenstate is partially separated from what we 
call an all- atom state by the fracture Jq S . The cut Jq S 
starts at 5 ~ — 17.8cj and K ~ 13.7 uj. For compar- 
ison, on the same figure we give the single GPE col- 
lapse boundary, j-^-y = — 0.67 u. We conclude that in 
the ground state the fracture lines correspond to the de- 
generacy points in the parameter space where two config- 
urations happen to have the same (lowest) total energy. 

In a simple AMBEC there is no formal collapse like 
the one that exists in a single GPE. The formal collapse 
is prevented by the two-atom resonance term in Eq. (16), 
which effectively counteracts the scattering term. In Ta- 
ble I we gave the fraction and the half-size of the atomic 
distribution along calculated numerically on both 
sides of the fracture for various values of the detun- 
ing. We saw a substantial decrease in the half-size of 
the atomic distribution across the fracture as one moves 
away from the origin, which may imitate a collapse as 
observed in a single GPE. In an AMBEC the collapse of 
the ground state becomes a controllable feature. We may 
create the collapse by driving the system across /gg, or 
we may drive the system around the fracture. The latter 
case should allow, we speculate, creation of a mixed or 
atom-molecule state of desired size in a continuous fash- 
ion and without any parametric excitations. 

B. Twin-State 

The twin state with its local eigenstates is shown in 
Fig. 13. The position of the fracture frw is deduced 
from /q S by extending the definition of Jq S to negative 
values of K. In the vicinity of Jtw (shaded region in 
Fig. 13) the mixed or atom-molecule state of the twin 
state disappears. As the nearby all-molecule eigenstate 
is unstable, the twin state decays rapidly after crossing 
the fracture. For that reason a conversion of atoms to 
molecules cannot take place in the twin state. In our 
numerical simulations we occasionally (not shown in this 
paper) observed a sequence of parametric excitations re- 
sembling a staircase in the shaded region. It appears as 
if the system went through a series of local eigenstates. 
However, we have not been able to develop a quantitative 
characterization of such states, if any. 

VII. DISCUSSION 

We have studied the eigenstates of the simple atom- 
molecule Bose-Einstein Condensate both numerically 
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FIG. 12: Map of the ground state with the fractures /q S and f% s separating the three participating eigenstates. 
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FIG. 13: Map of the twin state with the fracture frw separating the all-molecule state from the mixed or atom-molecule state. 
For non-infinitesimal fractions of the atoms, the onset of the instability occurs somewhere in the shaded region. 
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and analytically. Technically, we are investigating eigen- 
states of a nonlinear differential operator. The problem 
proves much trickier than analysis of eigenstates of lin- 
ear operators. To meet the challenge, we have devised 
the local eigenstate paradigm. The nature of a global 
eigenstate is qualitatively different in different regions 
of the problem parameters, and the transition from one 
type of local eigenstate to another is typically nonanalytic 
when a boundary between local eigenstates, a fracture, is 
crossed. When the problem parameters are varied in time 
across a fracture, the system will respond to the nonan- 
alyticity with breakdown of adiabaticity and parametric 
excitations. 

Two types of parametric excitations have been encoun- 
tered. First, the atomic condensate may start to oscillate 
erratically in shape while maintaining a more or less con- 
stant norm (ground state near /q S ). Second, the atomic 
condensate may degrade into seemingly chaotic behavior 
both in shape and in norm (twin state near fxw)- I n 
fact, a parameter sweep across the fracture frw is what 
is going on in the "Bosenova" experiments [30], and the 
ensuing loss of the atomic condensate makes a tempting 
parallel to the chaotic behavior we are predicting. As far 
as we can tell, in this particular experiment our picture 
of parametric excitations is as good as any theoretical 
description of the Bosenova offered to date [7, 31-33]. 

However, for reasons of time scales, the present ap- 
proach cannot be the panacea for all of the existing 
experiments. We are dealing with stationary or near- 
stationary states of atoms and molecules in a trap, so the 
trap frequency is a crucial parameter in our calculations. 
In contrast, say, the Ramsey fringe experiments in an 
atom-molecule condensate [30] probe time scales orders 
of magnitude shorter than the time of a back-and-forth 
oscillation of an atom in the trap. The trap frequency 
cannot be a significant factor in these experiments. Also, 
while the chaotic decay of the atomic condensate af- 
ter a sweep across Jtw is reminiscent of both modula- 
tional instability of an AMBEC [34] and the phenomenon 
we [35] have dubbed rogue dissociation, the discussions 
in Refs. [35-37] are for free atoms. The (likely nontrivial) 
connection between free and harmonically trapped gases 
remains to be worked out. 

Even though in the present experiments on AMBEC 
at least the atoms are usually trapped, often the trap 
is actually irrelevant in the dynamics. In contrast, we 
have argued that, in a trap, it might be possibly in prin- 
ciple to manipulate a condensate around the fracture 
Jq S , and thereby create a collapsed state of the atom- 
molecule condensate smoothly and in a completely con- 
trollable manner. Studies of a trapped AMBEC on the 
time scale of the trap, both experimental and theoretical, 
might have other such intriguing surprises in store. 
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APPENDIX A: NUMERICS 



All numerical calculations in this paper were performed 
using the novel DS-method [26] which we developed to in- 
tegrate systems of nonlinear Schrodinger equations such 
as Eqs. (3). We use pointwise representation for the am- 
plitudes and the operators, with N = 1024 being the 
total number of points and Rq = 11 the spatial extent 
of the system. The harmonic potential of the trap is 
isotropic with u) = 1. The trap Hamiltonians H a and 
H m contain only the s-wave parts, 



2 m a>m dr 



-4>, (Al) 



where m a — 1 and m m — 2 are the atomic and molecular 
masses. Integration is performed with a fixed time step 
chosen according to the prescription 



At 



0.01 



max 



1, y'maxt K(t), max t \S(t)\ 



(A2) 



In each step of the integration in complex time some loss 
of the norm of the atomic and molecular amplitudes oc- 
curs. The renormalization scheme wc used to compensate 
for the losses consists of the following formulae 



ip(ir + iAr) 
ijj{iT + iAr) 



l+n„ 



2n m +n a 



+n a 



■ ip(iT + iAr), 

■ tp(iT + iAr). 



(A3) 



Here n a = (<p(ir + iAr)\(p(iT + iAr)) and n m — 
(ip(iT + iAr)\ip{iT + iAr)) are the norms of the atomic 
and molecular amplitudes before renormalization. This 
renormalization furnishes a fixed point with tp and ip cor- 
rect to first order in the integration time step At. When 
integrating in real time, the amplitudes ip and ip may also 
gradually lose norm. In that case they are renormalized 
according to the same prescription, Eq. (A3). Half-size 
of the atomic distribution ip is a numerical solution R1/2 
of the following equation. 



f °° drr 2 \<p(r)\ 2 



(A4) 



We use Eq. (A4) to estimate the extent of the condensate 
in the trap, and in particular to study the oscillations the 
condensate may undergo following the parametric exci- 
tation at a fracture. 
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